Chapter 8　 Space Filling

8.1　 Introduction

This chapter focuses on the Space filling problem * 1 and explains Apollonian Gasket , which is one of the methods to solve it .

Since this chapter focuses on the algorithm explanation of Apollonius Gasket, it deviates a little from the story of graphic programming.

8.2　 Space filling problem

The space filling problem is the problem of finding a method to fill the inside of one closed plane as much as possible with a certain shape without overlapping. This problem is an area that has been studied for a long time, especially in the fields of geometry and combinatorial optimization. Since there are innumerable combinations of what kind of plane to fill with what kind of shape, various methods have been proposed for each combination.

To give a few examples

• Rectangular packing * 2　: O-Tree method
• Polygonal packing * 3 : Bottom-Left method
• Circle packing * 4　　: Apollonian Gasket
• Triangular packing * 5 : Sierpinski Gasket

And so on, there are many other techniques. In this chapter, we will explain about Apollonian Gasket among the above.

The Space filling problem is known to be NP-hard, and it is currently difficult to always fill the plane 100% with any of the above algorithms. The same is true for the Apollonian Gasket, which cannot completely fill the inside of a circle.

8.3　 Apollonian Gasket

The Apollonian Gasket is a type of fractal figure generated from three circles that touch each other. This is a kind of the earliest fractal figure, and it is said that one of the research results of plane geometry could be one solution of the Space filling problem, not the algorithm proposed to solve the Space filling problem. is. The name is named after Apollonius of Perga, a Greek mortar in BC.

First, assuming that the three circles that touch each other are C1, C2, and C3, respectively, Apollonius discovered that there are two non-intersecting circles C4 and C5 that touch all of C1, C2, and C3. These C4 and C5 are Apollonius circles for C1, C2 and C3 (details will be described later) .

Figure 8.1: C1, C2, C3 and the two circles C4, C5 in contact with it

Now, if we consider C4 as opposed to C1, C2, we can get two new Apollonius circles for C1, C2, and C4. Of these two circles, one will be C3 and the other will be the new circle C6.

Considering the Apollonius circles for all combinations, such as (C1, C2, C5), (C2, C3, C4), (C1, C3, C4), you can get at least one new circle for each. I can do it. By repeating this infinitely, a set of circles that touch each other is created. This set of circles is the Apollonian Gasket.

Figure 8.2: Apollonian Gasket

https://upload.wikimedia.org/wikipedia/commons/e/e6/Apollonian_gasket.svg

8.4　 Apollonian Gasket Calculation

From here, I will explain the calculation method of Apollonian Gasket in order while looking at the actual program. A sample program is available on Github, so please download it from there if necessary.

URL：https://github.com/IndieVisualLab/UnityGraphicsProgramming2

8.4.1　 Preparation

In programming the Apollonian Gasket, this time we have prepared our own class to represent the circle and a structure to handle complex numbers.

using UnityEngine;

public class Circle
{
    public float Curvature
    {
        get { return 1f / this.radius; }
    }
    public Complex Complex
    {
        get; private set;
    }
    public float Radius
    {
        get { return Mathf.Abs(this.radius); }
    }
    public Vector2 Position
    {
        get { return this.Complex.Vec2; }
    }

    private float radius = 0f;


    public Circle(Complex complex, float radius)
    {
        this.radius = radius;
        this.Complex = complex;
    }

    /// ...
    /// Below, a function to check the relationship between circles is implemented.
    /// Whether they are in contact, intersect, include, etc.
    /// ...
}

Part of the implementation of a class that represents a circle.

If you have basic programming knowledge, there should be nothing difficult. In addition, Complexis a complex number structure prepared by myself this time , and is Curvaturecalled curvature, both of which are necessary values ​​for calculating the Apollonian Gasket.

using UnityEngine;
using System;
using System.Globalization;

public struct Complex
{
    public static readonly Complex Zero = new Complex(0f, 0f);
    public static readonly Complex One = new Complex(1f, 0f);
    public static readonly Complex ImaginaryOne = new Complex(0f, 1f);

    public float Real
    {
        get { return this.real; }
    }
    public float Imaginary
    {
        get { return this.imaginary; }
    }
    public float Magnitude
    {
        get { return Abs(this); }
    }
    public float SqrMagnitude
    {
        get { return SqrAbs(this); }
    }
    public float Phase
    {
        get { return Mathf.Atan2(this.imaginary, this.real); }
    }
    public Vector2 Vec2
    {
        get { return new Vector2(this.real, this.imaginary); }
    }

    [SerializeField]
    private float real;
    [SerializeField]
    private float imaginary;


    public Complex(Vector2 vec2) : this(vec2.x, vec2.y) { }

    public Complex(Complex other) : this(other.real, other.imaginary) { }

    public Complex(float real, float imaginary)
    {
        this.real = real;
        this.imaginary = imaginary;
    }

    /// ...
    /// Below, the function to calculate the complex number is implemented.
    /// Four arithmetic operations, absolute value calculation, etc.
    /// ...
}

A structure for handling complex numbers.

C # has a Complexstructure, but it has been included since .Net 4.0. At the time of writing this chapter, Unity's .Net 4.6 support was in the Experimental stage, so I decided to prepare it myself.

8.4.2　 Calculation of the first three circles

As a prerequisite for calculating the Apollonian Gasket, there must be three circles tangent to each other. Therefore, in this program, three circles with randomly determined radii are generated, and the coordinates are calculated and arranged so that they touch each other.

private void CreateFirstCircles(
    out Circle c1, out Circle c2, out Circle c3)
{
    var r1 = Random.Range (
        this.firstRadiusMin, this.firstRadiusMax
    );
    var r2 = Random.Range (
        this.firstRadiusMin, this.firstRadiusMax
    );
    var r3 = Random.Range(
        this.firstRadiusMin, this.firstRadiusMax
    );

    // Get random coordinates
    var p1 = this.GetRandPosInCircle(
        this.fieldRadiusMin,
        this.fieldRadiusMax
    );
    c1 = new Circle(new Complex(p1), r1);

    // Calculate the center coordinates of the tangent circle based on p1
    var p2 = -p1.normalized * ((r1 - p1.magnitude) + r2);
    c2 = new Circle(new Complex(p2), r2);

    // Calculate the center coordinates of a circle tangent to two circles
    var p3 = this.GetThirdVertex(p1, p2, r1 + r2, r2 + r3, r1 + r3);
    c3 = new Circle(new Complex(p3), r3);
}

private Vector2 GetRandPosInCircle(float fieldMin, float fieldMax)
{
    // Get the right angle
    var theta = Random.Range (0f, Mathf.PI * 2f);

    // Calculate the appropriate distance
    var radius = Mathf.Sqrt(
        2f * Random.Range(
            0.5f * fieldMin * fieldMin,
            0.5f * fieldMax * fieldMax
        )
    );

    // Convert from polar coordinate system to Euclidean plane
    return new Vector2(
        radius * Mathf.Cos(theta),
        radius * Mathf.Sin(theta)
    );
}

private Vector2 GetThirdVertex(
    Vector2 p1, Vector2 p2, float rab, float rbc, float rca)
{
    var p21 = p2 - p1;

    // Calculate the angle by the cosine theorem
    var theta = Mathf.Acos(
        (rab * rab + rca * rca - rbc * rbc) / (2f * rca * rab)
    );

    // Calculate and add the angle of the starting point
    // theta is just an angle in the triangle, not an angle in the plane
    theta += Mathf.Atan2(p21.y, p21.x);

    // Add the coordinates converted from the polar coordinate system to the Euclidean plane to the starting coordinates.
    return p1 + new Vector2(
        rca * Mathf.Cos(theta),
        rca * Mathf.Sin(theta)
    );
}

CreateFirstCirclesBy calling the function, the initial condition of 3 yen is generated.

Randomly into three radial first r1,r2,r3decided, then GetRandPosInCirclethe function r1to determine the center coordinates of a circle having a radius (hereinafter C1). This function returns random coordinates inside a circle that is fieldMingreater fieldMaxthan or equal to the radius of the origin center .

Figure 8.3: Area where random coordinates are generated

Next, r2the center coordinates of a circle with a radius (below C2) to calculate. First , calculate the distance from the origin to the center of C2 by (r1-p1.magnitude) + r2 . By multiplying this by the normalized coordinates of C1 which is sign-inverted, r2the center coordinates of a circle with a radius adjacent to C1 can be obtained.

Figure 8.4: Position represented by the p2 vector

r3The center coordinates of the circle with the last radius (hereinafter C3) are GetThirdVertexcalculated by the function, but this calculation uses the cosine theorem . As most readers may have learned in high school, the cosine theorem is a theorem that holds between the length of the sides of a triangle and the cosine of an internal angle (cos). △ in ABC, a = BC, b = CA, c = AB, alpha = ∠CAB When

a^2 = c^2 + b^2 - 2cbcosα

Is the law of cosines.

Figure 8.5: Triangle ABC

You may be wondering why you need a triangle to think about the center of a circle, but in fact, the relationship between the three circles makes it possible to think of a very easy-to-use triangle. Considering a triangle whose apex is the center of C1, C2, and C3, since these three circles are in contact with each other, the length of each side of the triangle can be known from the radius of the circle.

Figure 8.6: Triangle ABC and Circles C1, C2, C3

When the law of cosines is transformed

cosα = \frac{c^2 + b^2 - a^2}{2cb}

Therefore, we can solve the cosine from the lengths of the three sides. Once the angle and distance between the two sides are found, the center coordinates of C3 can be found based on the center coordinates of C1. ..

You have now generated the three tangent circles you need as an initial condition.

8.4.3　 Calculation of the circle tangent to C1, C2, C3

Based on the three circles C1, C2, and C3 generated in the previous section, the circles tangent to them are calculated. Two parameters, radius and center coordinates, are required to create a new circle, so each is calculated.

radius

We first calculate from the radius, which can be calculated by Descartes's circle theorem . Descartes's circle theorem is that for four circles C1, C2, C3, C4 that touch each other, the curvature * 6 is k1, k2, k3, k4, respectively.

(k_1 + k_2 + k_3 + k_4) ^ 2 = 2 ({k_1} ^ 2 + {k_2} ^ 2 + {k_3} ^ 2 + {k_4} ^ 2)

Is true. This is a quadratic equation for the radii of four circles, but if you organize this equation

k_4 = k_1 + k_2 + k_3 \ pm 2 \ sqrt {k_1k_2 + k_2k_3 + k_3k_1}

If the three circles C1, C2, and C3 are known, the curvature of the fourth circle C4 can be obtained. Since the curvature is the reciprocal of the radius, the radius of the circle can be known by taking the reciprocal of the curvature.

Here, the curvature of C4 is obtained by compounding two, but one solution is always positive and the other is either positive or negative. When the curvature of C4 is positive, it is circumscribed to C1, C2, C3, and when it is negative, it is inscribed to C1, C2, C3 (including 3 circles). In other words, the fourth circle C4 can be considered in two patterns, and there is a possibility that both can be considered.

Figure 8.7: Positive and negative curvature

The following part is programming this.

// Curvature calculation
var k1 = this.Circle1.Curvature;
var k2 = this.Circle2.Curvature;
var k3 = this.Circle3.Curvature;

var plusK = k1 + k2 + k3 + 2f * Mathf.Sqrt (k1 * k2 + k2 * k3 + k3 * k1);
var minusK = k1 + k2 + k3 - 2f * Mathf.Sqrt (k1 * k2 + k2 * k3 + k3 * k1);

This Cartesian circle theorem was later rediscovered by a chemist named Sodi, and the C1, C2, C3, and C4 circles are called Sodi's circles.

Center coordinates

Next is the calculation of the center coordinates, which is calculated by the Cartesian complex number theorem , which has a shape similar to the Cartesian circle theorem . The Cartesian complex number theorem is that the center coordinates of the circles C1, C2, C3, C4 that touch each other on the complex plane are z1, z2, z3, z4, and the curvature is k1, k2, k3, k4.

(k_1z_1 + k_2z_2 + k_3z_3 + k_4z_4) ^ 2 = 2 ({k_1} ^ 2 {z_1} ^ 2 + {k_2} ^ 2 {z_2} ^ 2 + {z_3} ^ 2 {z_3} ^ 2 + {k_4} ^ 2 {z_4} ^ 2)

Is true. To organize this formula for z4

z4 = \frac{z_1k_1 + z_2k_2 + z_3k_3 \pm 2\sqrt{k_1k_2z_1z_2 + k_2k_3z_2z_3 + k_3k_1z_3z_1}}{k4}

Since it can be transformed into, the center coordinates of the circle C4 can be obtained with this.

Here, two curvatures were obtained when calculating the radius, but two can also be obtained by compounding the Cartesian complex number theorem. However, unlike the curvature calculation, one of the two is the correct Sodi circle, so you need to determine which is correct.

The following part is programming this.

/// Calculation of center coordinates
var ck1 = Complex.Multiply(this.Circle1.Complex, k1);
var ck2 = Complex.Multiply(this.Circle2.Complex, k2);
var ck3 = Complex.Multiply(this.Circle3.Complex, k3);

var plusZ = ck1 + ck2 + ck3
    + Complex.Multiply(Complex.Sqrt(ck1 * ck2 + ck2 * ck3 + ck3 * ck1), 2f);
var minusZ = ck1 + ck2 + ck3
    - Complex.Multiply(Complex.Sqrt(ck1 * ck2 + ck2 * ck3 + ck3 * ck1), 2f);

var recPlusK = 1f / plusK;
var recMinusK = 1f / minusK;

// Judgment of Sodi's circle
this.GetGasket(
    new Circle(Complex.Divide(plusZ, plusK), recPlusK),
    new Circle(Complex.Divide(minusZ, plusK), recPlusK),
    out c4
);

this.GetGasket(
    new Circle(Complex.Divide(plusZ, minusK), recMinusK),
    new Circle(Complex.Divide(minusZ, minusK), recMinusK),
    out c5
);

/// Judgment of Sodi's circle
(c1.IsCircumscribed(c4, CalculationAccuracy)
    || c1.IsInscribed(c4, CalculationAccuracy)) &&
(c2.IsCircumscribed(c4, CalculationAccuracy)
    || c2.IsInscribed(c4, CalculationAccuracy)) &&
(c3.IsCircumscribed(c4, CalculationAccuracy)
    || c3.IsInscribed(c4, CalculationAccuracy))

public bool IsCircumscribed(Circle c, float accuracy)
{
    var d = (this.Position - c.Position).sqrMagnitude;
    var abs = Mathf.Abs(d - Mathf.Pow(this.Radius + c.Radius, 2));

    return abs <= accuracy * accuracy;
}

public bool IsInscribed(Circle c, float accuracy)
{
    var d = (this.Position - c.Position).sqrMagnitude;
    var abs = Mathf.Abs(d - Mathf.Pow(this.Radius - c.Radius, 2));

    return abs <= accuracy * accuracy;
}

Now, based on the initial conditions C1, C2, C3, we have obtained two circles tangent to them (hereinafter C4, C5).

8.4.4　 Apollonian Gasket Calculation

At this point, you can easily calculate the Apollonian Gasket. Simply repeat the calculation performed in "8.4.3 Calculation of the circle tangent to C1, C2, C3" .

In the previous section, we found the circles C4 and C5 that are tangent to C1, C2 and C3 found in "8.4.2 Calculation of the first three circles" . Next, find the circles that touch (C1, C2, C4) (C1, C2, C5) (C2, C3, C4) (C2, C3, C5) (C3, C1, C4) (C3, C1, C5). I will continue.

Here, even if the circles are in contact with each other on the combination, they may actually overlap with other circles. Therefore, after determining whether it is the correct Sodi circle, it is also necessary to confirm that it does not overlap all the circles that have been requested so far.

Figure 8.8: Of the circles C7 and C8 tangent to C1, C4 and C6, C8 overlaps C2 and is not included in the Apollonian Gasket.

Then you will get a new circle that touches each one. After that, in the same way, we will continue to find new tangent circles for all combinations of the original circle to find the tangent circle and the newly found tangent circle.

Mathematically, the set of circles obtained by repeating this procedure infinitely is the Apollonian Gasket, but it is not possible to handle the program infinitely. Therefore, in this program, if the radius of the newly obtained tangent circle is less than a certain value, the condition that the processing is completed is given for the combination.

The following part is programming this.

private void Awake()
{
    // Generate the initial condition of 3 yen
    Circle c1, c2, c3;
    this.CreateFirstCircles(out c1, out c2, out c3);
    this.circles.Add(c1);
    this.circles.Add(c2);
    this.circles.Add(c3);

    this.soddys.Enqueue(new SoddyCircles(c1, c2, c3));

    while(this.soddys.Count > 0)
    {
        // Calculate Sodi's circle
        var soddy = this.soddys.Dequeue();

        Circle c4, c5;
        soddy.GetApollonianGaskets(out c4, out c5);

        this.AddCircle(c4, soddy);
        this.AddCircle(c5, soddy);
    }
}

private void AddCircle(Circle c, SoddyCircles soddy)
{
    if(c == null || c.Radius <= MinimumRadius)
    {
        return;
    }
    // If the curvature is negative, no questions asked and added
    // Circles with negative curvature appear only once
    else if(c.Curvature < 0f)
    {
        this.circles.Add(c);
        soddy.GetSoddyCircles(c).ForEach(s => this.soddys.Enqueue(s));

        return;
    }

    // Check if it covers other circles
    for(var i = 0; i < this.circles.Count; i++)
    {
        var o = this.circles[i];

        if(o.Curvature < 0f)
        {
            continue;
        }
        else if(o.IsMatch(c, CalculationAccuracy) == true)
        {
            return;
        }
    }

    this.circles.Add(c);
    soddy.GetSoddyCircles(c).ForEach(s => this.soddys.Enqueue(s));
}

You have now successfully requested the Apollonian Gasket.

Figure 8.9: Execution result on Unity

8.5　 Summary

So far, we have walked through the steps required to calculate the Apollonian Gasket. As I explained at the beginning, the Apollonian Gasket has a stronger meaning as a fractal figure.

However, if we remove the limitation of the plane this time and jump out into the world of space, it will become difficult to talk about it, and the meaning of filling (packing) from the fractal figure will become stronger. The proposition of sphere-packing space is a field that has been controversial for hundreds of years, including the existence of famous mathematical conjectures such as the Kepler conjecture.

The Space filling problem is also useful in practical terms. It is applied in a wide range of fields such as optimization of VLSI layout design, optimization of cutting out parts such as cloth, and automation and optimization of UV development.

This time, I chose the Apollonian Gasket, which is relatively easy to understand and interesting. If you are interested in packing itself, check out the algorithms introduced at the beginning.

Fill the inside of the object with the object. I think it can be used as a new expression method in unexpected places.

8.6　 Reference

• https://ja.wikipedia.org/wiki/ Apollonian Gasket
• https://ja.wikipedia.org/wiki/ Descartes' Circle Theorem
• http://paulbourke.net/fractals/randomtile/